1 Abstract

The goal of this study is to look at territory size of the Louisiana Waterthrush during breeding season at three locations; Wall Branch stream located in Rotary Park in Clarksville Tennessee 37040, Dry Creek located in Cheatham Wildlife Management Area in Ashland City Tennessee 37015, Big Hollow located in Beaman Park in Ashland City Tennessee 37015. Mapping these territories in these streams will show not only the territories of the Louisiana Waterthrush present but also possible areas of concentrated use that may aid in the location of their nesting sites.

2 Feild Methods

We first mapped out our study area of the streams (a few weeks before the arrival of the first birds) with a handheld global positioning system (GPS) so we had a reliable map of our study sites when out in the field, the streams chosen for this project were not well mapped before and Google Maps did not provide a detailed enough map that would be useful. To map the stream, we used a handheld GPS unit the Garmin Etrex-30 model to mark waypoints every 3-5 meters apart. While this was done, the way point number was written down as well stream information pertaining to that way-point location. Stream information included; stream width, stream type (confluence or tributary), and noted if any water was present. The notation of water is important because of the feeding habits of the Louisiana Waterthrush, having this information about suitable habitats will be helpful later when searching out individuals. Maps were then made using ArcMap. Sampling for territory mapping occurred within the first week of the Louisiana Waterthrush being spot at the stream locations, streams were checked in rotation every three days until a Waterthrush was spotted. Sampling period was from March 21, 2019 until May 10, 2019. This period of time worked well because we were able to get location points before birds had paired up as well as when the nests were established, and the first clutches of the season were laid. Location points had to be collected in a way that wouldn’t affect the focal individual, but also allowed time for the observer to observe all birds in a location. Observations were made primarily in the morning, since this is when the birds sing most often, although there were a few afternoon observations due to scheduling interferences. It was important that the observer recording location points did not scare the bird during this process. The Louisiana Waterthrush spends a large majority of its time foraging in the stream, so location points where taken via direct observation. The points were taken once the Louisiana Waterthrush had left the area and was a safe distance away that taking the point would not scare the bird away. This meant that the observer needed to make a mental note of any points needed to be taken until the bird was more than 15m away. Location points were taken where the bird was first encounter, then every 5 minutes within a 20-minute period, points were also taken where the bird sang, when the bird moved at least five meters as well as where the bird turned around in the stream (possibly marking the end of a territory). This usually meant following the bird while they foraged on the stream. The 20-minute period was chosen to make sure that there was still enough time to observe other birds in the section of the stream because of time constraints during the weekdays. Even though there way a 20-minute observational period, not all focal individuals were observed for the full 20 minutes; in that case the individual was considered “lost” (see protocol below). Occasionally extra location points were taken on individuals that were observed for longer outside the 20-minutes period, if time allowed it.

3 Kernal Density Maps

We will be making KDE maps so we will need to register with google maps to use their satelite version of our map area. The link below is where a lot of the following code came from. https://cfss.uchicago.edu/notes/raster-maps-with-ggmap/ You must register for google maps to get an AI (I used the one year free trial for this and they don’t auto charge once the year is over). OpenMapSource is another option as well if your java supports this. Once you have registered create your base map using google maps address or coordinates.

3.1 Basemap for Dry Creek,TN

Next, we can put our LOWA encounter cooridinates on the basemap we just created.

DCmap<-ggmap(basemap)+ geom_point(data=DryCreek, aes(Longitude, Latitude)) +
  geom_point(data = DryCreek, aes(Longitude, Latitude), size = 0.1) +
  theme(axis.title = element_text(face="bold")) + labs(x="Longitude", y="Latitude") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme_bw()
plot(DCmap)

Now it’s time to create KDE and MCP, but first we have to format a few things. It will be important to first put your data into SpatialPoints we will be using CRS (coordinate reference system) which is the format for most GPS data. The link below is where a lot of the following code came from. https://mhallwor.github.io/_pages/activities_GenerateTerritories

DCpts <- sp::SpatialPoints(coords = cbind(DryCreek$Longitude, DryCreek$Latitude),
                                proj4string = sp::CRS("+init=epsg:4326"))
head(DCpts)

## class       : SpatialPoints 
## features    : 1 
## extent      : -87.18876, -87.18876, 36.18876, 36.18876  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

## put all the data into a single SpatialPointsDataFrame (spdf)
DC_spdf <- sp::SpatialPointsDataFrame(DCpts, DryCreek)
head(DC_spdf)

#although we don't have seperate species we need to let R know that all points are for one stream system and this will create a list of one
DC_sep <- split(x = DC_spdf, f = DC_spdf$Name, drop = FALSE)

Time to make the MCP.

DCmcp <- lapply(DC_sep, FUN = function(x){rgeos::gConvexHull(x)})

##this makes polygon from the list of one when we told R that all variables are from one stream system
DCmcp <- mapply(DCmcp, names(DCmcp), 
                  SIMPLIFY = FALSE,
                  FUN = function(x,y){x@polygons[[1]]@ID <- y
                  return(x)})

DCmcp <- do.call(rbind,DCmcp)
DCmcp <- SpatialPolygonsDataFrame(Sr = DCmcp,
                                   data = data.frame(Bird = names(DCmcp)),
                                   match.ID = FALSE)
plot(DCmcp)

Now that we have an MCP for the stream, onto KDE.

## Step one: do least squares cross-validation to estimate bandwidth (you may get a warning message but keep going)
bw <- lapply(DC_sep, FUN = function(x){ks::Hlscv(x@coords)})

## Warning in ks::Hlscv(x@coords): Data contain duplicated values: LSCV is not
## well-behaved in this case

## Step two: generate kde

DC_kde <-mapply(DC_sep,bw,
                     SIMPLIFY = FALSE,
                     FUN = function(x,y){
                       raster(kde(x@coords,h=y))})
# This code makes a custom function called getContour. 
# Inputs:
#    kde = kernel density estimate
#    prob = probabily - default is 0.95

getContour <- function(kde, prob = 0.95){
  # set all values 0 to NA
  kde[kde == 0]<-NA
  # create a vector of raster values
  kde_values <- raster::getValues(kde)
  # sort values 
  sortedValues <- sort(kde_values[!is.na(kde_values)],decreasing = TRUE)
  # find cumulative sum up to ith location
  sums <- cumsum(as.numeric(sortedValues))
  # binary response is value in the probabily zone or not
  p <- sum(sums <= prob * sums[length(sums)])
  # Set values in raster to 1 or 0
  kdeprob <- raster::setValues(kde, kde_values >= sortedValues[p])
  # return new kde
  return(kdeprob)}

DC_95kde <- lapply(DC_kde,
                    FUN = getContour,prob = 0.95)

These next plots put MCP on top of KDE

plot(DC_kde[[1]])+
plot(DCmcp[1,],add = TRUE)

## integer(0)

Time to overlap KDE onto the base map. Code found from http://data-analytics.net/cep/Schedule_files/geospatial.html

KDEmap<-DCmap +
  stat_density2d(aes(x = Longitude, y = Latitude, fill = ..level..,alpha=..level..), bins = 10, geom = "polygon", data = DryCreek) +
  scale_fill_gradient(low = "red", high = "green")+
  ggtitle("Dry Creek, TN")

3.1.1 Final KDE Map

Compleated KDE Heatmap

4 KDE Map for Big Hollow, TN

Big_Hollow <- read_csv("/Users/Kirstenglish/Desktop/Big Hollow.csv")

## Parsed with column specification:
## cols(
##   Name = col_character(),
##   Latitude = col_double(),
##   Longitude = col_double(),
##   Number = col_double()
## )

View(Big_Hollow)

4.1 Basemap

basemap<- get_map(location = "Ridgetop Trail
Ashland City, TN 37015", zoom=15)

## Source : https://maps.googleapis.com/maps/api/staticmap?center=Ridgetop%20Trail%0AAshland%20City,%20TN%2037015&zoom=15&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Ridgetop+Trail%0AAshland+City,+TN+37015&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg

ggmap(basemap)

Put our LOWA encounter cooridinates on the basemap we just created.

#time to put points on the map
BHmap<-ggmap(basemap)+ geom_point(data=Big_Hollow, aes(Longitude, Latitude)) +
  geom_point(data = Big_Hollow, aes(Longitude, Latitude), size = 0.1) +
  theme(axis.title = element_text(face="bold")) + labs(x="Longitude", y="Latitude") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme_bw()
plot(BHmap)

Put data into spatial points like before to be able to create the KDE and MCP. we will be using CRS (coordinate reference system) which is the format for most GPS data.

Time to make the MCP.

BHmcp <- lapply(BH_sep, FUN = function(x){rgeos::gConvexHull(x)})

##this makes polygon from the list of one when we told R that all variables are from one stream system
BHmcp <- mapply(BHmcp, names(BHmcp), 
                SIMPLIFY = FALSE,
                FUN = function(x,y){x@polygons[[1]]@ID <- y
                return(x)})

BHmcp <- do.call(rbind,BHmcp)
BHmcp <- SpatialPolygonsDataFrame(Sr = BHmcp,
                                  data = data.frame(Bird = names(BHmcp)),
                                  match.ID = FALSE)
plot(BHmcp)

Now that we have an MCP for the stream, onto KDE.

bw <- lapply(BH_sep, FUN = function(x){ks::Hlscv(x@coords)})
## Step two: generate kde

BH_kde <-mapply(BH_sep,bw,
                SIMPLIFY = FALSE,
                FUN = function(x,y){
                  raster(kde(x@coords,h=y))})

# This code makes a custom function called getContour. 
# Inputs:
#    kde = kernel density estimate
#    prob = probabily - default is 0.95

getContour <- function(kde, prob = 0.95){
  # set all values 0 to NA
  kde[kde == 0]<-NA
  # create a vector of raster values
  kde_values <- raster::getValues(kde)
  # sort values 
  sortedValues <- sort(kde_values[!is.na(kde_values)],decreasing = TRUE)
  # find cumulative sum up to ith location
  sums <- cumsum(as.numeric(sortedValues))
  # binary response is value in the probabily zone or not
  p <- sum(sums <= prob * sums[length(sums)])
  # Set values in raster to 1 or 0
  kdeprob <- raster::setValues(kde, kde_values >= sortedValues[p])
  # return new kde
  return(kdeprob)}

BH_95kde <- lapply(BH_kde,
                   FUN = getContour,prob = 0.95)

These next plots put MCP on top of KDE

plot(BH_kde[[1]])+
  plot(BHmcp[1,],add = TRUE)

## integer(0)

Time to overlap KDE onto the base map. Code found from http://data-analytics.net/cep/Schedule_files/geospatial.html

4.1.1 Final KDE Map

Compleated KDE Heatmap

5 Kernal Density Map for Wall Branch, TN

Create the basemap for the site location. This site is Wall Branch stream in Tennesse

basemap<- get_map((location = " 2561 Alex Overlook Way
Clarksville, TN 37043, 36.497818, -87.267428"), zoom=15, maptype = "satellite")

## Source : https://maps.googleapis.com/maps/api/staticmap?center=%202561%20Alex%20Overlook%20Way%0AClarksville,%20TN%2037043,%2036.497818,%20-87.267428&zoom=15&size=640x640&scale=2&maptype=satellite&language=en-EN&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg

## Source : https://maps.googleapis.com/maps/api/geocode/json?address=2561+Alex+Overlook+Way%0AClarksville,+TN+37043,+36.497818,+-87.267428&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg

ggmap(basemap)

LOWA encounter cooridinates on the basemap we just created.

#time to put points on the map
WBmap<-ggmap(basemap)+ geom_point(data=Wall_Branch, aes(Longitude, Latitude)) +
  geom_point(data = Wall_Branch, aes(Longitude, Latitude), color= "black", size = 0.1) +
  theme(axis.title = element_text(face="bold")) + labs(x="Longitude", y="Latitude") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme_classic()
plot(WBmap)

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_point).

WBpts <- sp::SpatialPoints(coords = cbind(Wall_Branch$Longitude, Wall_Branch$Latitude),
                           proj4string = sp::CRS("+init=epsg:4326"))
head(WBpts)

## class       : SpatialPoints 
## features    : 1 
## extent      : -87.26923, -87.26923, 36.49896, 36.49896  (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

WB_spdf <- sp::SpatialPointsDataFrame(WBpts, Wall_Branch)
head(WB_spdf)

WB_sep <- split(x = WB_spdf, f = WB_spdf$Name, drop = FALSE)

Time to make the MCP.

WBmcp <- lapply(WB_sep, FUN = function(x){rgeos::gConvexHull(x)})

##this makes polygon from the list of one when we told R that all variables are from one stream system
WBmcp <- mapply(WBmcp, names(WBmcp), 
                SIMPLIFY = FALSE,
                FUN = function(x,y){x@polygons[[1]]@ID <- y
                return(x)})
WBmcp <- do.call(rbind,WBmcp)
WBmcp <- SpatialPolygonsDataFrame(Sr = WBmcp,
                                  data = data.frame(Bird = names(WBmcp)),
                                  match.ID = FALSE)
plot(WBmcp)

Now that we have an MCP for the stream, onto KDE.

bw <- lapply(WB_sep, FUN = function(x){ks::Hlscv(x@coords)})

## Step two: generate kde
WB_kde <-mapply(WB_sep,bw,
                SIMPLIFY = FALSE,
                FUN = function(x,y){
                  raster(kde(x@coords,h=y))})
# This code makes a custom function called getContour. 
# Inputs:
#    kde = kernel density estimate
#    prob = probabily - default is 0.95

getContour <- function(kde, prob = 0.95){
  # set all values 0 to NA
  kde[kde == 0]<-NA
  # create a vector of raster values
  kde_values <- raster::getValues(kde)
  # sort values 
  sortedValues <- sort(kde_values[!is.na(kde_values)],decreasing = TRUE)
  # find cumulative sum up to ith location
  sums <- cumsum(as.numeric(sortedValues))
  # binary response is value in the probabily zone or not
  p <- sum(sums <= prob * sums[length(sums)])
  # Set values in raster to 1 or 0
  kdeprob <- raster::setValues(kde, kde_values >= sortedValues[p])
  # return new kde
  return(kdeprob)}

WB_95kde <- lapply(WB_kde,
                   FUN = getContour,prob = 0.95)

These next plots put MCP on top of KDE

## These next plots put MCP on top of KDE, (territory for this map is very linear)
plot(WB_kde[[1]])+
  plot(WBmcp[1,],add = TRUE)

## integer(0)

Time to overlap KDE onto the base map. Code found from http://data-analytics.net/cep/Schedule_files/geospatial.html

WBkde<-WBmap +
  stat_density2d(aes(x = Longitude, y = Latitude, fill = ..level..,alpha=..level..), bins = 10, geom = "polygon", data = Wall_Branch) +
  scale_fill_gradient(low = "red", high = "green")+
  ggtitle("Wall Branch, TN")

Compleated KDE Heatmap

6 Species Distrubution Map of the Louisiana Waterthrush

library(dismo)

## 
## Attaching package: 'dismo'

## The following object is masked from 'package:ggmap':
## 
##     geocode

library(rgbif)
library(rdryad)

## Registered S3 method overwritten by 'crul':
##   method                 from
##   as.character.form_file httr

library(readxl)
library(ggridges)

## 
## Attaching package: 'ggridges'

## The following object is masked from 'package:ggplot2':
## 
##     scale_discrete_manual

library(viridis)

## Loading required package: viridisLite

library(rasterVis)

## Loading required package: lattice

## Loading required package: latticeExtra

## Loading required package: RColorBrewer

## 
## Attaching package: 'latticeExtra'

## The following object is masked from 'package:ggplot2':
## 
##     layer

extent <- extent(-102,-80,10,45)

LOWA_dismo <- gbif("seiurus", species = "motacilla", ext = extent,
             geo = TRUE, sp = TRUE, download = TRUE,
             removeZeros = TRUE)

## 942 records found

## 0-300-600-900-942 records downloaded

LOWA_xy <- as.data.frame(cbind(LOWA_dismo@coords[,1],LOWA_dismo@coords[,2]))
colnames(LOWA_xy) <- c("longitude","latitude")

us <- map_data("world")
ggplot(data = LOWA_xy, aes(x=longitude, y=latitude)) +
  geom_polygon(data = us, aes(x=long, y = lat, group = group),
               fill = "white", color="black") +
  geom_point() + xlab("Longitude") + ylab("Latitude") +
  coord_fixed(xlim = c(-102,-80), ylim = c(10,45))

bioclim <- getData(name = "worldclim", res = 2.5, var = "bio")

names(bioclim) <- c("Ann Mean Temp","Mean Diurnal Range","Isothermality","Temperature Seasonality","Max Temp Warmest Mo","Min Temp Coldest Mo","Ann Temp Range","Mean Temp Wettest Qtr","Mean Temp Driest Qtr","Mean Temp Warmest Qtr","Mean Temp Coldest Qtr","Annual Precip","Precip Wettest Mo","Precip Driest Mo","Precip Seasonality","Precip Wettest Qtr","Precip Driest Qtr","Precip Warmest Qtr","Precip Coldest Qtr")

bio_extent <- extent(x = c(
  min(LOWA_xy$longitude),
  max(LOWA_xy$longitude),
  min(LOWA_xy$latitude),
  max(LOWA_xy$latitude)))

bioclim_extent <- crop(x = bioclim, y = bio_extent)
bioclim_model <- bioclim(x = bioclim_extent, p = LOWA_xy)
presence_model <- dismo::predict(object = bioclim_model, 
                                 x = bioclim_extent, 
                                 ext = bio_extent)

gplot(presence_model) + 
  geom_raster(aes(fill=value)) +
  geom_polygon(data = us, aes(x= long, y = lat, group = group),
               fill = NA, color="black") +
  geom_point(data = LOWA_xy, aes(x = longitude, y = latitude), color = "black", alpha = 0.5) +
  scale_fill_gradientn(colours=c("brown","yellow","darkgreen"), "Probability") +
  coord_fixed(xlim = c(-102,-80), ylim = c(10,45)) +
  xlab("Longitude") + ylab("Latitude") + ggtitle("Probability of LOWA Occurrence") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + theme(legend.position = "right") +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

7 Feather Under Confocal Microscope

Compleated KDE Heatmap